Baby Mental Life: Study 3
Preregistered analyses
2018-05-07
# load required libraries
library(tidyverse)
library(langcog) # source: https://github.com/langcog/langcog-package
library(psych)
library(lme4)
library(kableExtra)
library(ggdendro)
library(dendextend)
# set theme for ggplots
theme_set(theme_bw())
chosen_rot <- "oblimin"# run source code (extra home-made functions)
source("./scripts/max_factors_efa.R")
source("./scripts/plot_fun.R")
source("./scripts/reten_fun.R")
source("./scripts/data_prep.R")NAs introduced by coercionattributes are not identical across measure variables;
they will be droppedJoining, by = "question_qualtrics"
Joining, by = "question"
Joining, by = c("ResponseId", "Q36", "Q36_coded_KH", "Q36_coded_LSK", "Q36_coded_KW", "Q36_coded_KH_simple", "Q36_coded_LSK_simple", "Q36_coded_KW_simple", "MATCH_simple")
Joining, by = c("ResponseId", "attn_free_coded")
Column `ResponseId` joining character vector and factor, coercing into character vectorColumn `attn_free_coded` joining character vector and factor, coercing into character vectorJoining, by = "ResponseId"
Joining, by = "ResponseId"
NAs introduced by coercionJoining, by = "capacity"
Column `capacity` joining character vector and factor, coercing into character vectorJoining, by = "capacity"
Column `capacity` joining character vector and factor, coercing into character vectorJoining, by = "capacity"
Column `capacity` joining character vector and factor, coercing into character vectorJoining, by = "capacity"
Column `capacity` joining character vector and factor, coercing into character vector“Baby Mental Life: Study 3” was conducted on MTurk on 2019-04-24.
Our planned sample was 300 participants; based on Study 2, we initially recruited 352 participants. After filtering out participants who failed at least one of our attention checks, we ended up retaining 279 participants (retention rate: 85.5%). We then supplemented by recruiting an additional 28 participants, with the goal of attaining our planned sample (300 - 279 = 21 participants, 279/352 * 28 = 22). At each stage, we recruited women and men through separate studies, in hopes of acquiring a roughly equal split between genders.
In the end, we ended up with a sample of 301 participants who passed our attention checks, 252 of whom came from unique GPS coordinates (83.7%).
For this first pass, these data include participants where there is another participant with an identical set of GPS coordinates as recorded by Qualtrics.
Each participant assessed children’s mental capacities at 13 target ages between the ages of 0 and 5 years. For each target, they rated 8 mental capacities on a scale from 0 (not at all capable) to 100 (completely capable). In contrast to Study 2, participants completed all 13 ratings of a particular mental capacity on a single “trial” (rather than completing all 8 ratings of a particular target age on a single trial). In addition, participants answered questions about the “developmental factors” that might contribute to development in this domain.
For more details about the study, see our preregistration here.
Study 1 EFA
# load in EFA results from study 1
efa_S1 <- readRDS("../study 1/s1_efa.rds")heatmap_fun(efa_S1) +
labs(title = paste0("STUDY 1 Parallel Analysis (rotation: ", chosen_rot, ")"),
subtitle = "'% var.' indicates the amount of shared variance explained (total = 100%)")Joining, by = "capacity"
Joining, by = "factor"
Plots: Capacity ratings as a function of target age and domain
ggplot(d_cap_rating %>%
arrange(ResponseId, domain, capacity, target_ord) %>%
mutate(domain = recode_factor(
domain,
"BOD" = "Bodily sensations",
"NEG" = "Negative emotions",
"POS" = "Social abilities & positive emotions",
"COG" = "Cognition & control",
.default = NA_character_),
capacity = gsub("_", " ", capacity)),
aes(x = target_num, y = response, group = ResponseId,
color = domain)) +
facet_wrap(~ domain ~ capacity, ncol = 4) +
geom_path(alpha = 0.1) +
geom_smooth(aes(group = capacity),
method = "loess",
color = "black") +
scale_x_continuous(breaks = as.numeric(levels(factor(d_cap_rating$target_num))),
labels = levels(d_cap_rating$target_ord)) +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
labs(title = "Perceptions of development, by domain",
subtitle = "Black line = loess smoothing",
x = "Target age (numeric)",
y = "Response (0 = not at all capable, 100 = completely capable)")
ggplot(d_cap_rating %>%
arrange(ResponseId, domain, capacity, target_ord) %>%
mutate(domain = recode_factor(
domain,
"BOD" = "Bodily sensations",
"NEG" = "Negative emotions",
"POS" = "Social abilities & positive emotions",
"COG" = "Cognition & control",
.default = NA_character_),
capacity = gsub("_", " ", capacity)),
aes(x = sqrt(target_num), y = response, group = ResponseId,
color = domain)) +
facet_wrap(~ domain ~ capacity, ncol = 4) +
geom_path(alpha = 0.1) +
geom_smooth(aes(group = capacity),
method = "loess",
color = "black") +
scale_x_continuous(breaks = sqrt(as.numeric(levels(factor(d_cap_rating$target_num)))),
labels = levels(d_cap_rating$target_ord)) +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
labs(title = "Perceptions of development, by domain",
subtitle = "Black line = loess smoothing",
x = "Target age (numeric, square-root transformed)",
y = "Response (0 = not at all capable, 100 = completely capable)")
ggplot(d_cap_rating %>%
arrange(ResponseId, domain, capacity, target_ord) %>%
mutate(domain = recode_factor(
domain,
"BOD" = "Bodily sensations",
"NEG" = "Negative emotions",
"POS" = "Social abilities & positive emotions",
"COG" = "Cognition & control",
.default = NA_character_),
capacity = gsub("_", " ", capacity)),
aes(x = target_ord, y = response, group = ResponseId,
color = domain)) +
facet_wrap(~ domain ~ capacity, ncol = 4) +
geom_path(alpha = 0.1) +
geom_smooth(aes(group = capacity),
method = "loess",
color = "black") +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
labs(title = "Perceptions of development, by domain",
subtitle = "Black line = loess smoothing",
x = "Target age (ordinal)",
y = "Response (0 = not at all capable, 100 = completely capable)")
Regression models, version 1: Capacity ratings as a function of target age and domain
contrasts_domain_eff_noCOG <- cbind(BOD = c(1, -1, 0, 0),
NEG = c(0, -1, 1, 0),
POS = c(0, -1, 0, 1))
contrasts_domain_eff_noPOS <- cbind(BOD = c(1, 0, 0, -1),
NEG = c(0, 0, 1, -1),
COG = c(0, 1, 0, -1))
contrasts_domain_eff_noNEG <- cbind(BOD = c(1, 0, -1, 0),
POS = c(0, 0, -1, 1),
COG = c(0, 1, -1, 0))
contrasts_domain_eff_noBOD <- cbind(POS = c(-1, 0, 0, 1),
NEG = c(-1, 0, 1, 0),
COG = c(-1, 1, 0, 0))Target age in months
Effect-coded (comparing all domains to the grand mean, except “bodily sensations”)
Linear effects only
# r1_noBOD <- lmer(response ~ target_num * domain
# + (target_num + domain | ResponseId)
# + (target_num | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12),
# contrasts = list(domain = contrasts_domain_eff_noBOD))
# saveRDS(r1_noBOD, "./models/r1_noBOD.RDS")
r1_noBOD <- readRDS("./models/r1_noBOD.RDS")
summary(r1_noBOD, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula:
response ~ target_num * domain + (target_num + domain | ResponseId) +
(target_num | capacity)
Data: d_cap_rating %>% mutate(target_num = target_num/12)
REML criterion at convergence: 275913.1
Scaled residuals:
Min 1Q Median 3Q Max
-5.2212 -0.4875 0.0406 0.5273 4.6405
Random effects:
Groups Name Variance Std.Dev. Corr
ResponseId (Intercept) 115.676 10.755
target_num 9.926 3.151 -0.49
domainCOG 403.013 20.075 -0.55 -0.02
domainNEG 402.083 20.052 0.19 -0.49 0.28
domainPOS 260.980 16.155 -0.05 -0.41 0.52 0.51
capacity (Intercept) 203.077 14.251
target_num 7.030 2.651 -0.89
Residual 348.192 18.660
Number of obs: 31304, groups: ResponseId, 301; capacity, 8
Fixed effects:
Estimate Std. Error t value
(Intercept) 55.0228 5.1044 10.780
target_num 8.1378 0.9572 8.502
domainPOS -3.5835 8.7464 -0.410
domainNEG 9.0183 8.7652 1.029
domainCOG -44.0243 8.7650 -5.023
target_num:domainPOS 3.7284 1.6278 2.291
target_num:domainNEG -2.2729 1.6278 -1.396
target_num:domainCOG 5.4199 1.6278 3.330Non-linear effects
# r2_noBOD <- lmer(response ~ poly(target_num, 3) * domain
# + (target_num + domain | ResponseId)
# + (target_num | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12),
# contrasts = list(domain = contrasts_domain_eff_noBOD))
# saveRDS(r2_noBOD, "./models/r2_noBOD.RDS")
r2_noBOD <- readRDS("./models/r2_noBOD.RDS")
summary(r2_noBOD, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(target_num, 3) * domain + (target_num + domain |
ResponseId) + (target_num | capacity)
Data: d_cap_rating %>% mutate(target_num = target_num/12)
REML criterion at convergence: 271075.5
Scaled residuals:
Min 1Q Median 3Q Max
-5.2790 -0.4882 0.0098 0.5135 5.3651
Random effects:
Groups Name Variance Std.Dev. Corr
ResponseId (Intercept) 118.031 10.864
target_num 10.123 3.182 -0.49
domainCOG 406.943 20.173 -0.55 -0.02
domainNEG 406.006 20.150 0.18 -0.48 0.28
domainPOS 264.904 16.276 -0.06 -0.40 0.52 0.51
capacity (Intercept) 203.382 14.261
target_num 7.044 2.654 -0.89
Residual 297.181 17.239
Number of obs: 31304, groups: ResponseId, 301; capacity, 8
Fixed effects:
Estimate Std. Error t value
(Intercept) 66.506 3.973 16.739
poly(target_num, 3)1 2270.326 267.286 8.494
poly(target_num, 3)2 -929.790 17.239 -53.936
poly(target_num, 3)3 445.691 17.239 25.854
domainPOS 1.678 6.808 0.246
domainNEG 5.811 6.832 0.851
domainCOG -36.376 6.831 -5.325
poly(target_num, 3)1:domainPOS 1040.177 454.392 2.289
poly(target_num, 3)2:domainPOS -937.701 29.859 -31.405
poly(target_num, 3)3:domainPOS 559.497 29.859 18.738
poly(target_num, 3)1:domainNEG -634.099 454.392 -1.395
poly(target_num, 3)2:domainNEG 165.960 29.859 5.558
poly(target_num, 3)3:domainNEG -87.169 29.859 -2.919
poly(target_num, 3)1:domainCOG 1512.077 454.392 3.328
poly(target_num, 3)2:domainCOG 74.583 29.859 2.498
poly(target_num, 3)3:domainCOG -189.126 29.859 -6.334
Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need itEffect-coded (comparing all domains to the grand mean, except “negative emotions”)
Linear effects only
# r1_noNEG <- lmer(response ~ target_num * domain
# + (target_num + domain | ResponseId)
# + (target_num | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12),
# contrasts = list(domain = contrasts_domain_eff_noNEG))
# saveRDS(r1_noNEG, "./models/r1_noNEG.RDS")
r1_noNEG <- readRDS("./models/r1_noNEG.RDS")
summary(r1_noNEG, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula:
response ~ target_num * domain + (target_num + domain | ResponseId) +
(target_num | capacity)
Data: d_cap_rating %>% mutate(target_num = target_num/12)
REML criterion at convergence: 275913.1
Scaled residuals:
Min 1Q Median 3Q Max
-5.2212 -0.4875 0.0406 0.5273 4.6405
Random effects:
Groups Name Variance Std.Dev. Corr
ResponseId (Intercept) 115.676 10.755
target_num 9.926 3.151 -0.49
domainCOG 403.022 20.075 -0.55 -0.02
domainNEG 402.081 20.052 0.19 -0.49 0.28
domainPOS 260.982 16.155 -0.05 -0.41 0.52 0.51
capacity (Intercept) 203.390 14.261
target_num 7.039 2.653 -0.89
Residual 348.192 18.660
Number of obs: 31304, groups: ResponseId, 301; capacity, 8
Fixed effects:
Estimate Std. Error t value
(Intercept) 55.0228 5.1082 10.772
target_num 8.1378 0.9578 8.497
domainBOD 38.5896 8.7600 4.405
domainPOS -3.5835 8.7531 -0.409
domainCOG -44.0243 8.7717 -5.019
target_num:domainBOD -6.8755 1.6288 -4.221
target_num:domainPOS 3.7284 1.6288 2.289
target_num:domainCOG 5.4199 1.6288 3.328Non-linear effects
# r2_noNEG <- lmer(response ~ poly(target_num, 3) * domain
# + (target_num + domain | ResponseId)
# + (target_num | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12),
# contrasts = list(domain = contrasts_domain_eff_noNEG))
# saveRDS(r2_noNEG, "./models/r2_noNEG.RDS")
r2_noNEG <- readRDS("./models/r2_noNEG.RDS")
summary(r2_noNEG, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(target_num, 3) * domain + (target_num + domain |
ResponseId) + (target_num | capacity)
Data: d_cap_rating %>% mutate(target_num = target_num/12)
REML criterion at convergence: 271075.5
Scaled residuals:
Min 1Q Median 3Q Max
-5.2790 -0.4882 0.0098 0.5135 5.3651
Random effects:
Groups Name Variance Std.Dev. Corr
ResponseId (Intercept) 118.030 10.864
target_num 10.124 3.182 -0.49
domainCOG 406.935 20.173 -0.55 -0.02
domainNEG 406.006 20.150 0.18 -0.48 0.28
domainPOS 264.902 16.276 -0.06 -0.40 0.52 0.51
capacity (Intercept) 203.225 14.256
target_num 7.039 2.653 -0.89
Residual 297.181 17.239
Number of obs: 31304, groups: ResponseId, 301; capacity, 8
Fixed effects:
Estimate Std. Error t value
(Intercept) 66.506 3.972 16.745
poly(target_num, 3)1 2270.326 267.205 8.497
poly(target_num, 3)2 -929.790 17.239 -53.936
poly(target_num, 3)3 445.691 17.239 25.854
domainBOD 28.888 6.814 4.240
domainPOS 1.678 6.805 0.247
domainCOG -36.376 6.829 -5.327
poly(target_num, 3)1:domainBOD -1918.156 454.249 -4.223
poly(target_num, 3)2:domainBOD 697.158 29.859 23.349
poly(target_num, 3)3:domainBOD -283.202 29.859 -9.485
poly(target_num, 3)1:domainPOS 1040.177 454.249 2.290
poly(target_num, 3)2:domainPOS -937.701 29.859 -31.405
poly(target_num, 3)3:domainPOS 559.497 29.859 18.738
poly(target_num, 3)1:domainCOG 1512.077 454.249 3.329
poly(target_num, 3)2:domainCOG 74.583 29.859 2.498
poly(target_num, 3)3:domainCOG -189.126 29.859 -6.334
Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need itEffect-coded (comparing all domains to the grand mean, except “cognition & control”)
Linear effects only
# r1_noCOG <- lmer(response ~ target_num * domain
# + (target_num + domain | ResponseId)
# + (target_num | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12),
# contrasts = list(domain = contrasts_domain_eff_noCOG))
# saveRDS(r1_noCOG, "./models/r1_noCOG.RDS")
r1_noCOG <- readRDS("./models/r1_noCOG.RDS")
summary(r1_noCOG, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula:
response ~ target_num * domain + (target_num + domain | ResponseId) +
(target_num | capacity)
Data: d_cap_rating %>% mutate(target_num = target_num/12)
REML criterion at convergence: 275913.1
Scaled residuals:
Min 1Q Median 3Q Max
-5.2212 -0.4875 0.0406 0.5273 4.6405
Random effects:
Groups Name Variance Std.Dev. Corr
ResponseId (Intercept) 115.677 10.755
target_num 9.926 3.151 -0.49
domainCOG 403.014 20.075 -0.55 -0.02
domainNEG 402.089 20.052 0.19 -0.49 0.28
domainPOS 260.979 16.155 -0.05 -0.41 0.52 0.51
capacity (Intercept) 203.084 14.251
target_num 7.030 2.651 -0.89
Residual 348.192 18.660
Number of obs: 31304, groups: ResponseId, 301; capacity, 8
Fixed effects:
Estimate Std. Error t value
(Intercept) 55.0228 5.1044 10.779
target_num 8.1378 0.9572 8.502
domainBOD 38.5896 8.7534 4.409
domainNEG 9.0183 8.7653 1.029
domainPOS -3.5835 8.7466 -0.410
target_num:domainBOD -6.8755 1.6278 -4.224
target_num:domainNEG -2.2729 1.6278 -1.396
target_num:domainPOS 3.7284 1.6278 2.290Non-linear effects
# r2_noCOG <- lmer(response ~ poly(target_num, 3) * domain
# + (target_num + domain | ResponseId)
# + (target_num | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12),
# contrasts = list(domain = contrasts_domain_eff_noCOG))
# saveRDS(r2_noCOG, "./models/r2_noCOG.RDS")
r2_noCOG <- readRDS("./models/r2_noCOG.RDS")
summary(r2_noCOG, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(target_num, 3) * domain + (target_num + domain |
ResponseId) + (target_num | capacity)
Data: d_cap_rating %>% mutate(target_num = target_num/12)
REML criterion at convergence: 271075.5
Scaled residuals:
Min 1Q Median 3Q Max
-5.2790 -0.4882 0.0098 0.5135 5.3651
Random effects:
Groups Name Variance Std.Dev. Corr
ResponseId (Intercept) 118.029 10.864
target_num 10.124 3.182 -0.49
domainCOG 406.933 20.173 -0.55 -0.02
domainNEG 406.008 20.150 0.18 -0.48 0.28
domainPOS 264.906 16.276 -0.06 -0.40 0.52 0.51
capacity (Intercept) 203.498 14.265
target_num 7.046 2.654 -0.89
Residual 297.181 17.239
Number of obs: 31304, groups: ResponseId, 301; capacity, 8
Fixed effects:
Estimate Std. Error t value
(Intercept) 66.506 3.974 16.734
poly(target_num, 3)1 2270.326 267.336 8.492
poly(target_num, 3)2 -929.790 17.239 -53.936
poly(target_num, 3)3 445.691 17.239 25.854
domainBOD 28.888 6.818 4.237
domainNEG 5.811 6.834 0.850
domainPOS 1.678 6.809 0.246
poly(target_num, 3)1:domainBOD -1918.156 454.480 -4.221
poly(target_num, 3)2:domainBOD 697.158 29.859 23.349
poly(target_num, 3)3:domainBOD -283.202 29.859 -9.485
poly(target_num, 3)1:domainNEG -634.099 454.480 -1.395
poly(target_num, 3)2:domainNEG 165.960 29.859 5.558
poly(target_num, 3)3:domainNEG -87.169 29.859 -2.919
poly(target_num, 3)1:domainPOS 1040.177 454.480 2.289
poly(target_num, 3)2:domainPOS -937.701 29.859 -31.405
poly(target_num, 3)3:domainPOS 559.497 29.859 18.738
Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need itTarget age in months (square-root transformation)
NOTE: All of these models fail to meet our criterion of random effects being correlation < 0.90. Because of this, I have not included in this print-out and I have not yet modeled non-linear effects after square-root transformation. The next step will be to try simpler random effects models (with fewer random slopes).
Regression models, version 2: Capacity ratings as a function of target age, for each domain separately
NOTE: All of the preregistered models fail to meet our criterion of random effects being correlation < 0.90. Because of this, I instead used simpler random effects models (with fewer random slopes), as indicated in the preregistration.
Target age in months
Bodily sensations
Linear effects only
# r5a_BOD <- lmer(response ~ target_num
# + (1 | ResponseId)
# + (1 | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12) %>%
# filter(domain == "BOD"))
# saveRDS(r5a_BOD, "./models/r5a_BOD.RDS")
r5a_BOD <- readRDS("./models/r5a_BOD.RDS")
summary(r5a_BOD, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula: response ~ target_num + (1 | ResponseId) + (1 | capacity)
Data:
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==
"BOD")
REML criterion at convergence: 59314.9
Scaled residuals:
Min 1Q Median 3Q Max
-8.2674 -0.1305 0.0997 0.2482 3.8964
Random effects:
Groups Name Variance Std.Dev.
ResponseId (Intercept) 97.626 9.881
capacity (Intercept) 1.135 1.066
Residual 101.015 10.051
Number of obs: 7826, groups: ResponseId, 301; capacity, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 93.61241 0.95668 97.85
target_num 1.26233 0.07205 17.52Non-linear effects
# r6a_BOD <- lmer(response ~ poly(target_num, 3)
# + (1 | ResponseId)
# + (1 | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12) %>%
# filter(domain == "BOD"))
# saveRDS(r6a_BOD, "./models/r6a_BOD.RDS")
r6a_BOD <- readRDS("./models/r6a_BOD.RDS")
summary(r6a_BOD, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula:
response ~ poly(target_num, 3) + (1 | ResponseId) + (1 | capacity)
Data:
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==
"BOD")
REML criterion at convergence: 59090.3
Scaled residuals:
Min 1Q Median 3Q Max
-8.1734 -0.2154 0.0249 0.3262 4.1567
Random effects:
Groups Name Variance Std.Dev.
ResponseId (Intercept) 97.728 9.886
capacity (Intercept) 1.136 1.066
Residual 98.365 9.918
Number of obs: 7826, groups: ResponseId, 301; capacity, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 95.3937 0.9514 100.262
poly(target_num, 3)1 176.0851 9.9179 17.754
poly(target_num, 3)2 -116.3162 9.9179 -11.728
poly(target_num, 3)3 81.2443 9.9179 8.192Negative emotions
Linear effects only
# r5a_NEG <- lmer(response ~ target_num
# + (1 | ResponseId)
# + (1 | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12) %>%
# filter(domain == "NEG"))
# saveRDS(r5a_NEG, "./models/r5a_NEG.RDS")
r5a_NEG <- readRDS("./models/r5a_NEG.RDS")
summary(r5a_NEG, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula: response ~ target_num + (1 | ResponseId) + (1 | capacity)
Data:
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==
"NEG")
REML criterion at convergence: 72902.7
Scaled residuals:
Min 1Q Median 3Q Max
-3.8109 -0.5686 0.0777 0.7164 3.1800
Random effects:
Groups Name Variance Std.Dev.
ResponseId (Intercept) 476.5 21.83
capacity (Intercept) 241.9 15.55
Residual 576.7 24.01
Number of obs: 7826, groups: ResponseId, 301; capacity, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 64.0411 11.0765 5.782
target_num 5.8649 0.1721 34.069Non-linear effects
# r6a_NEG <- lmer(response ~ poly(target_num, 3)
# + (1 | ResponseId)
# + (1 | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12) %>%
# filter(domain == "NEG"))
# saveRDS(r6a_NEG, "./models/r6a_NEG.RDS")
r6a_NEG <- readRDS("./models/r6a_NEG.RDS")
summary(r6a_NEG, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula:
response ~ poly(target_num, 3) + (1 | ResponseId) + (1 | capacity)
Data:
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==
"NEG")
REML criterion at convergence: 72561.3
Scaled residuals:
Min 1Q Median 3Q Max
-3.6464 -0.6310 0.0716 0.6993 3.4952
Random effects:
Groups Name Variance Std.Dev.
ResponseId (Intercept) 477.4 21.85
capacity (Intercept) 241.9 15.55
Residual 553.1 23.52
Number of obs: 7826, groups: ResponseId, 301; capacity, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 72.32 11.07 6.530
poly(target_num, 3)1 818.11 23.52 34.785
poly(target_num, 3)2 -381.91 23.52 -16.239
poly(target_num, 3)3 179.26 23.52 7.622Negative emotions
Linear effects only
# r5a_COG <- lmer(response ~ target_num
# + (1 | ResponseId)
# + (1 | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12) %>%
# filter(domain == "COG"))
# saveRDS(r5a_COG, "./models/r5a_COG.RDS")
r5a_COG <- readRDS("./models/r5a_COG.RDS")
summary(r5a_COG, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula: response ~ target_num + (1 | ResponseId) + (1 | capacity)
Data:
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==
"COG")
REML criterion at convergence: 67265.1
Scaled residuals:
Min 1Q Median 3Q Max
-5.1158 -0.6063 -0.0316 0.5361 4.0924
Random effects:
Groups Name Variance Std.Dev.
ResponseId (Intercept) 254.66 15.958
capacity (Intercept) 22.37 4.729
Residual 279.59 16.721
Number of obs: 7826, groups: ResponseId, 301; capacity, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 10.9985 3.4776 3.163
target_num 13.5577 0.1199 113.104Non-linear effects
# r6a_COG <- lmer(response ~ poly(target_num, 3)
# + (1 | ResponseId)
# + (1 | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12) %>%
# filter(domain == "COG"))
# saveRDS(r6a_COG, "./models/r6a_COG.RDS")
r6a_COG <- readRDS("./models/r6a_COG.RDS")
summary(r6a_COG, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula:
response ~ poly(target_num, 3) + (1 | ResponseId) + (1 | capacity)
Data:
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==
"COG")
REML criterion at convergence: 66491.6
Scaled residuals:
Min 1Q Median 3Q Max
-4.8966 -0.6171 -0.0328 0.5891 4.1837
Random effects:
Groups Name Variance Std.Dev.
ResponseId (Intercept) 255.68 15.99
capacity (Intercept) 22.37 4.73
Residual 253.16 15.91
Number of obs: 7826, groups: ResponseId, 301; capacity, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 30.130 3.474 8.673
poly(target_num, 3)1 1891.202 15.911 118.860
poly(target_num, 3)2 -427.604 15.911 -26.874
poly(target_num, 3)3 128.282 15.911 8.062Target age in months (square-root transformation)
Bodily sensations
Linear effects only
# r7a_BOD <- lmer(response ~ sqrt(target_num)
# + (1 | ResponseId)
# + (1 | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12) %>%
# filter(domain == "BOD"))
# saveRDS(r7a_BOD, "./models/r7a_BOD.RDS")
r7a_BOD <- readRDS("./models/r7a_BOD.RDS")
summary(r7a_BOD, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula: response ~ sqrt(target_num) + (1 | ResponseId) + (1 | capacity)
Data:
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==
"BOD")
REML criterion at convergence: 59164.5
Scaled residuals:
Min 1Q Median 3Q Max
-8.1883 -0.2061 0.0803 0.2887 4.0990
Random effects:
Groups Name Variance Std.Dev.
ResponseId (Intercept) 97.702 9.884
capacity (Intercept) 1.136 1.066
Residual 99.036 9.952
Number of obs: 7826, groups: ResponseId, 301; capacity, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 92.0086 0.9643 95.42
sqrt(target_num) 3.5037 0.1628 21.53Non-linear effects
# r8a_BOD <- lmer(response ~ poly(sqrt(target_num), 3)
# + (1 | ResponseId)
# + (1 | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12) %>%
# filter(domain == "BOD"))
# saveRDS(r8a_BOD, "./models/r8a_BOD.RDS")
r8a_BOD <- readRDS("./models/r8a_BOD.RDS")
summary(r8a_BOD, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(sqrt(target_num), 3) + (1 | ResponseId) + (1 |
capacity)
Data:
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==
"BOD")
REML criterion at convergence: 59012.2
Scaled residuals:
Min 1Q Median 3Q Max
-8.0365 -0.2069 -0.0311 0.3508 4.3998
Random effects:
Groups Name Variance Std.Dev.
ResponseId (Intercept) 97.767 9.888
capacity (Intercept) 1.136 1.066
Residual 97.350 9.867
Number of obs: 7826, groups: ResponseId, 301; capacity, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 95.3937 0.9515 100.255
poly(sqrt(target_num), 3)1 214.2262 9.8666 21.712
poly(sqrt(target_num), 3)2 -106.5941 9.8666 -10.804
poly(sqrt(target_num), 3)3 38.9803 9.8666 3.951Negative emotions
Linear effects only
# r7a_NEG <- lmer(response ~ sqrt(target_num)
# + (1 | ResponseId)
# + (1 | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12) %>%
# filter(domain == "NEG"))
# saveRDS(r7a_NEG, "./models/r7a_NEG.RDS")
r7a_NEG <- readRDS("./models/r7a_NEG.RDS")
summary(r7a_NEG, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula: response ~ sqrt(target_num) + (1 | ResponseId) + (1 | capacity)
Data:
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==
"NEG")
REML criterion at convergence: 72630.4
Scaled residuals:
Min 1Q Median 3Q Max
-3.7869 -0.6042 0.0780 0.6935 3.4966
Random effects:
Groups Name Variance Std.Dev.
ResponseId (Intercept) 477.3 21.85
capacity (Intercept) 241.9 15.55
Residual 556.3 23.59
Number of obs: 7826, groups: ResponseId, 301; capacity, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 57.9857 11.0802 5.233
sqrt(target_num) 14.8338 0.3857 38.455Non-linear effects
# r8a_NEG <- lmer(response ~ poly(sqrt(target_num), 3)
# + (1 | ResponseId)
# + (1 | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12) %>%
# filter(domain == "NEG"))
# saveRDS(r8a_NEG, "./models/r8a_NEG.RDS")
r8a_NEG <- readRDS("./models/r8a_NEG.RDS")
summary(r8a_NEG, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(sqrt(target_num), 3) + (1 | ResponseId) + (1 |
capacity)
Data:
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==
"NEG")
REML criterion at convergence: 72538.6
Scaled residuals:
Min 1Q Median 3Q Max
-3.6105 -0.6399 0.0577 0.6911 3.6371
Random effects:
Groups Name Variance Std.Dev.
ResponseId (Intercept) 477.5 21.85
capacity (Intercept) 241.9 15.55
Residual 551.5 23.48
Number of obs: 7826, groups: ResponseId, 301; capacity, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 72.32 11.07 6.530
poly(sqrt(target_num), 3)1 906.98 23.48 38.622
poly(sqrt(target_num), 3)2 -190.43 23.48 -8.109
poly(sqrt(target_num), 3)3 -31.35 23.48 -1.335Cognition & control
Linear effects only
# r7a_COG <- lmer(response ~ sqrt(target_num)
# + (1 | ResponseId)
# + (1 | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12) %>%
# filter(domain == "COG"))
# saveRDS(r7a_COG, "./models/r7a_COG.RDS")
r7a_COG <- readRDS("./models/r7a_COG.RDS")
summary(r7a_COG, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula: response ~ sqrt(target_num) + (1 | ResponseId) + (1 | capacity)
Data:
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==
"COG")
REML criterion at convergence: 66686.1
Scaled residuals:
Min 1Q Median 3Q Max
-4.7865 -0.6315 -0.0466 0.5966 3.9492
Random effects:
Groups Name Variance Std.Dev.
ResponseId (Intercept) 255.45 15.98
capacity (Intercept) 22.37 4.73
Residual 258.94 16.09
Number of obs: 7826, groups: ResponseId, 301; capacity, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.3958 3.4832 -0.114
sqrt(target_num) 31.5957 0.2632 120.054Non-linear effects
# r8a_COG <- lmer(response ~ poly(sqrt(target_num), 3)
# + (1 | ResponseId)
# + (1 | capacity),
# d_cap_rating %>%
# mutate(target_num = target_num/12) %>%
# filter(domain == "COG"))
# saveRDS(r8a_COG, "./models/r8a_COG.RDS")
r8a_COG <- readRDS("./models/r8a_COG.RDS")
summary(r8a_COG, corr = F)Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(sqrt(target_num), 3) + (1 | ResponseId) + (1 |
capacity)
Data:
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==
"COG")
REML criterion at convergence: 66507.4
Scaled residuals:
Min 1Q Median 3Q Max
-4.8256 -0.6153 -0.0351 0.5816 4.1614
Random effects:
Groups Name Variance Std.Dev.
ResponseId (Intercept) 255.66 15.99
capacity (Intercept) 22.37 4.73
Residual 253.70 15.93
Number of obs: 7826, groups: ResponseId, 301; capacity, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 30.130 3.474 8.673
poly(sqrt(target_num), 3)1 1931.843 15.928 121.286
poly(sqrt(target_num), 3)2 128.315 15.928 8.056
poly(sqrt(target_num), 3)3 -153.107 15.928 -9.612Plots: Ratings for developmental factors, by domain and capacity
ggplot(d_dev_factor_rating %>%
mutate(domain = recode_factor(
domain,
"BOD" = "Bodily sensations",
"NEG" = "Negative emotions",
"POS" = "Social abilities & positive emotions",
"COG" = "Cognition & control",
.default = NA_character_)) %>%
mutate_at(vars(capacity, dev_factor),
funs(gsub("_", " ", .))) %>%
mutate(dev_factor = factor(
dev_factor,
levels = gsub("_", " ", levels(d_dev_factor_rating$dev_factor)))),
aes(x = dev_factor, y = response, color = domain)) +
facet_wrap(~ domain ~ capacity, ncol = 4) +
geom_jitter(alpha = 0.025, height = 0.2, width = 0.2) +
geom_pointrange(data = . %>% group_by(domain, capacity, dev_factor) %>%
multi_boot_standard(col = "response", na.rm = T),
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
color = "black", fatten = 4) +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
labs(title = "Ratings of developmental factors (Version 1)",
subtitle = "Error bars are bootstrapped 95% CIs",
x = "Developmental factor",
y = "Response (0 = plays no role, 6 = plays a very important role)")
ggplot(d_dev_factor_rating %>%
mutate(domain = recode_factor(
domain,
"BOD" = "Bodily sensations",
"NEG" = "Negative emotions",
"POS" = "Social abilities & positive emotions",
"COG" = "Cognition & control",
.default = NA_character_)) %>%
mutate_at(vars(capacity, dev_factor),
funs(gsub("_", " ", .))) %>%
mutate(dev_factor = factor(
dev_factor,
levels = gsub("_", " ", levels(d_dev_factor_rating$dev_factor)))),
aes(x = reorder(capacity, as.numeric(domain)),
y = response, color = domain)) +
facet_wrap(~ dev_factor, ncol = 5) +
geom_jitter(alpha = 0.025, height = 0.2, width = 0.2) +
geom_pointrange(data = . %>% group_by(domain, capacity, dev_factor) %>%
multi_boot_standard(col = "response", na.rm = T),
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
color = "black", fatten = 4) +
theme(legend.position = "top",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
guides(color = guide_legend(override.aes = list(alpha = 1, size = 3))) +
labs(title = "Ratings of developmental factors (Version 2)",
subtitle = "Error bars are bootstrapped 95% CIs",
color = "Domain",
x = "Capacity (by domain)",
y = "Response (0 = plays no role, 6 = plays a very important role)")
ggplot(d_dev_factor_rating %>%
mutate(dev_factor = factor(
gsub("_", " ", dev_factor),
levels = gsub("_", " ", levels(d_dev_factor_rating$dev_factor)))),
aes(x = dev_factor, y = response)) +
geom_jitter(alpha = 0.01, height = 0.3, width = 0.3, color = "blue") +
geom_pointrange(data = . %>% group_by(dev_factor) %>%
multi_boot_standard(col = "response", na.rm = T),
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
color = "black", fatten = 2) +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
labs(title = "Ratings of developmental factors (Version 3)",
subtitle = "Error bars are bootstrapped 95% CIs",
x = "Developmental factor",
y = "Response (0 = plays no role, 6 = plays a very important role)")
Exploration: EFA/clustering of developmental factors
It could be nice to reduce the number of “devleopmental factors” from 10 down to something more manageable for analysis. Here I explore EFA and hierarchical clustering as dimension reduction methods, and re-plot with these results in mind.
Parallel to what we’ve done for dimensionality reduction across capacities, to do this I’ll ignore some aspects of the within-subjects design of this study. I’ll consider each combination of participant and capacity to be a separate set of observations for this analysis (effectively, acting as though each participant were really 8 different participants, each of whom evaluated a different capacity on the 10 developmental factors included in the study).
d_dev_factor_rating_wide <- d_dev_factor_rating %>%
unite(key, ResponseId, domain, capacity) %>%
# mutate(response = response - 3) %>% # center at 0
spread(dev_factor, response) %>%
column_to_rownames("key")Exploratory factor analysis
reten_fun(d_dev_factor_rating_wide, chosen_rot)[1] 2fa.parallel(d_dev_factor_rating_wide)Parallel analysis suggests that the number of factors = 4 and the number of components = 2 
VSS(d_dev_factor_rating_wide, rotate = chosen_rot)convergence not obtained in GPFoblq. 1000 iterations used.convergence not obtained in GPFoblq. 1000 iterations used.
Very Simple Structure
Call: vss(x = x, n = n, rotate = rotate, diagonal = diagonal, fm = fm,
n.obs = n.obs, plot = plot, title = title, use = use, cor = cor)
VSS complexity 1 achieves a maximimum of 0.82 with 1 factors
VSS complexity 2 achieves a maximimum of 0.9 with 2 factors
The Velicer MAP achieves a minimum of 0.05 with 2 factors
BIC achieves a minimum of NA with 5 factors
Sample Size adjusted BIC achieves a minimum of NA with 5 factors
Statistics by number of factors 
Different factor retention protocols suggest retaining different numbers of factors (2, 4, or 5):
dev_factor_rating_efa2 <- fa(d_dev_factor_rating_wide,
nfactors = 2,
rotate = chosen_rot)heatmap_fun(dev_factor_rating_efa2) +
guides(fill = guide_colorbar(barheight = 8, barwidth = 0.5)) +
labs(title = "EFA of developmental factor ratings",
subtitle = "Retaining 2 factors (retention protocol: Weisman et al., 2017)")Joining, by = "capacity"
Joining, by = "factor"
dev_factor_rating_efa4 <- fa(d_dev_factor_rating_wide,
nfactors = 4,
rotate = chosen_rot)heatmap_fun(dev_factor_rating_efa4) +
guides(fill = guide_colorbar(barheight = 8, barwidth = 0.5)) +
labs(title = "EFA of developmental factor ratings",
subtitle = "Retaining 4 factors (retention protocol: parallel analysis)")Joining, by = "capacity"
Joining, by = "factor"
dev_factor_rating_efa5 <- fa(d_dev_factor_rating_wide,
nfactors = 5,
rotate = chosen_rot)heatmap_fun(dev_factor_rating_efa5) +
guides(fill = guide_colorbar(barheight = 8, barwidth = 0.5)) +
labs(title = "EFA of developmental factor ratings",
subtitle = "Retaining 5 factors (retention protocol: minimizing BIC)")Joining, by = "capacity"
Joining, by = "factor"
The 2- and 4-factor solutions each seem to offer some advantages… the 5-factor solution is sensible too, it just doesn’t offer much in the way of reducing the dataset!
Hierarchical clustering
dev_factor_rating_hclust <- d_dev_factor_rating_wide %>%
t() %>%
dist() %>%
hclust()
# dev_factor_rating_hclust_order <- data.frame(
# order = as.numeric(dev_factor_rating_hclust$order),
# dev_factor = as.character(dev_factor_rating_hclust$labels))
# cannot figure out how to do this automatically!
dev_factor_rating_hclust_order <- c("experiments", "people_teach",
"brain_changes", "observes_objects",
"observes_people", "interacts_people",
"preprogrammed", "womb_experiences",
"body_grows", "senses_improve")dev_factor_rating_hclust %>%
ggdendrogram(rotate = F) +
theme_minimal() +
theme(#axis.title = element_blank(),
# axis.text.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()) +
labs(title = "Hierarchical agglomerative clustering",
subtitle = "Complete linkage (default for stats::hclust() function)",
x = "Developmental factor", y = "Height")
# as.dendrogram() %>%
# set("labels_col", k = 2) %>%
# plot(horiz = T, xlim = c(-100, 200))This looks a lot like the 2-factor EFA solution to me!
Re-plotting by cluster
Let’s proceed with the clustering of developmental factors suggested by the 2-factor EFA solution and by hierarchical clustering.
d_dev_factor_rating_byclust <- d_dev_factor_rating %>%
mutate(dev_factor_cluster = case_when(
dev_factor %in% c("experiments", "people_teach", "brain_changes",
"observes_objects", "observes_people",
"interacts_people") ~ "extrinsic",
dev_factor %in% c("preprogrammed", "womb_experiences", "body_grows",
"senses_improve") ~ "intrinsic",
TRUE ~ NA_character_)) %>%
mutate(dev_factor_cluster = factor(dev_factor_cluster),
dev_factor = factor(dev_factor,
levels = dev_factor_rating_hclust_order),
response_cent = response - 3)
contrasts_cluster_eff <- cbind(EXT = c(1, -1))ggplot(d_dev_factor_rating_byclust %>%
mutate(domain = recode_factor(
domain,
"BOD" = "Bodily sensations",
"NEG" = "Negative emotions",
"POS" = "Social abilities & positive emotions",
"COG" = "Cognition & control",
.default = NA_character_)) %>%
mutate_at(vars(capacity, dev_factor),
funs(gsub("_", " ", .))) %>%
mutate(dev_factor = factor(
dev_factor,
levels = gsub("_", " ", dev_factor_rating_hclust_order))),
aes(x = dev_factor,
y = response, color = domain, shape = dev_factor_cluster)) +
facet_wrap(~ domain ~ capacity, ncol = 4) +
geom_jitter(alpha = 0.025, height = 0.2, width = 0.2, show.legend = F) +
geom_pointrange(data = . %>% group_by(domain, capacity,
dev_factor_cluster, dev_factor) %>%
multi_boot_standard(col = "response", na.rm = T),
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
color = "black", fatten = 4) +
theme(legend.position = "top",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
labs(title = "Ratings of developmental factors (Version 1)",
subtitle = "Error bars are bootstrapped 95% CIs",
shape = "Type of developmental factor (per EFA/clustering)",
x = "Developmental factor",
y = "Response (0 = plays no role, 6 = plays a very important role)")
ggplot(d_dev_factor_rating_byclust %>%
mutate(domain = recode_factor(
domain,
"BOD" = "Bodily sensations",
"NEG" = "Negative emotions",
"POS" = "Social abilities & positive emotions",
"COG" = "Cognition & control",
.default = NA_character_)) %>%
mutate_at(vars(capacity),
funs(gsub("_", " ", .))),
aes(x = reorder(capacity, as.numeric(domain)),
y = response, color = domain)) +
facet_wrap(dev_factor_cluster ~ dev_factor, ncol = 6) +
geom_jitter(alpha = 0.025, height = 0.2, width = 0.2) +
geom_pointrange(data = . %>% group_by(domain, capacity,
dev_factor_cluster, dev_factor) %>%
multi_boot_standard(col = "response", na.rm = T),
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
color = "black", fatten = 4) +
theme(legend.position = "top",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
guides(color = guide_legend(override.aes = list(alpha = 1, size = 3))) +
labs(title = "Ratings of developmental factors (Version 3)",
subtitle = "Error bars are bootstrapped 95% CIs",
color = "Domain",
x = "Capacity (by domain)",
y = "Response (0 = plays no role, 6 = plays a very important role)")
r_dev_factor_rating <- lmer(response ~ domain * dev_factor_cluster +
# (domain + dev_factor_cluster | ResponseId) +
# NOTE: including both random intercepts
# yields corr > 0.9
(1 + dev_factor_cluster | ResponseId) +
(1 | capacity) + (1 | dev_factor),
d_dev_factor_rating_byclust,
contrasts = list(
domain = contrasts_domain_eff_noBOD,
dev_factor_cluster = contrasts_cluster_eff))summary(r_dev_factor_rating)Linear mixed model fit by REML ['lmerMod']
Formula:
response ~ domain * dev_factor_cluster + (1 + dev_factor_cluster |
ResponseId) + (1 | capacity) + (1 | dev_factor)
Data: d_dev_factor_rating_byclust
REML criterion at convergence: 94588.4
Scaled residuals:
Min 1Q Median 3Q Max
-4.2100 -0.7160 0.0053 0.7145 3.4784
Random effects:
Groups Name Variance Std.Dev. Corr
ResponseId (Intercept) 0.78302 0.8849
dev_factor_clusterintrinsic 0.33241 0.5765 -0.19
dev_factor (Intercept) 0.46271 0.6802
capacity (Intercept) 0.09071 0.3012
Residual 2.80683 1.6754
Number of obs: 24080, groups:
ResponseId, 301; dev_factor, 10; capacity, 8
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.15782 0.24942 12.661
domainPOS 0.56248 0.18542 3.034
domainNEG -0.22905 0.18542 -1.235
domainCOG 0.21295 0.18542 1.148
dev_factor_clusterEXT 0.34626 0.22045 1.571
domainPOS:dev_factor_clusterEXT 0.33709 0.01909 17.662
domainNEG:dev_factor_clusterEXT -0.11557 0.01909 -6.055
domainCOG:dev_factor_clusterEXT 0.85384 0.01909 44.737
Correlation of Fixed Effects:
(Intr) dmnPOS dmnNEG dmnCOG d__EXT dPOS:_ dNEG:_
domainPOS 0.000
domainNEG 0.000 -0.333
domainCOG 0.000 -0.333 -0.333
dv_fctr_EXT -0.178 0.000 0.000 0.000
dmPOS:__EXT 0.000 -0.021 0.007 0.007 0.000
dmNEG:__EXT 0.000 0.007 -0.021 0.007 0.000 -0.333
dmCOG:__EXT 0.000 0.007 0.007 -0.021 0.000 -0.333 -0.333(Note that domain was effect-coded and leaves out the bodily sensations domain here, but we could re-code it in all the different ways to get the contrast between bodily sensations and the grand mean as desired.)
I see three take-aways from these regression results, all of which I think I can see reflected in the plots above:
- Participants perceived the domain of social connection and positive emotions to be particularly dramatically shaped by the set of developmental factors we included in this study (relative to the grand mean collapsing across domains).
- On the whole, participants endorsed “extrinsic” and “intrinsic” developmental factors equally strongly. If we were mapping this onto some sort of nature-nuruture distinction (not perfect, but tempting), I think we’d say that there’s no evidence for participants favoring either nature or nurture over the other as an explanation for growth in the capacities we included in this study.
- However, their endorsement of “extrinsic” and “intrinsic” developmental factors varied by domain. They perceived “extrinsic” factors to play a larger role in the domains of social connection and positive emotions and especially cognition and control, and a relatively small role in the domain of negative emotions. (Looking at the plots, it seems like extrinsic factors were perceived as much more important than intrinsic factors for cognition and control, a bit more important than intrinsic factors for social connection and positive emotions, roughly equally important to intrinsic factors for negative emotions, and less important than intrinsic factors for bodily sensations. These interaction terms just confirm that the difference between “extrinsic” and “intrinsic” varied significantly across domains.)
Table: Answers to (optional) free-response question about other developmental factors
d_dev_factor_other %>%
mutate(domain = recode_factor(
domain,
"BOD" = "Bodily sensations",
"NEG" = "Negative emotions",
"POS" = "Social abilities & positive emotions",
"COG" = "Cognition & control",
.default = NA_character_),
capacity = gsub("_", " ", capacity)) %>%
filter(response != "idk") %>%
distinct(domain, capacity, response) %>%
arrange(domain, capacity, response) %>%
kable(caption = "Free responses to optional question about other developmental factors that might play a role (excluding answers like 'N/A,' 'no,' 'not that I can think of', etc.)") %>%
kable_styling() %>%
collapse_rows(1:2)| domain | capacity | response |
|---|---|---|
| Bodily sensations | feeling pain | Again, you're born with this ability |
| all children develop differently so i rate as ive experienced | ||
| I believe that most living things can feel pain | ||
| I think experiencing physical pain is possibly as soon as a brain develops, long before birth | ||
| I think feeling pain is innate | ||
| I think pain is one of the first sensations that humans are able to actually comprehend | ||
| It is just a natural sensation. Our body is programmed to feel bad to avoid further injury or notifiy us that one has occured | ||
| Nervous system | ||
| None come to mind | ||
| Not that I can think of. I think everyone is capable of feeling pain, even at birth, as it is something physical. Whether they understand why they feel it though is a different story | ||
| Not that I know of | ||
| Pain is likely the very first thing we all experience when we are born. Childbirth is not exactly pleasant for the mother or the baby | ||
| Social interaction | ||
| The pain response is pre-programmed but, I feel it may take a couple of days at least for it to recognize it is pain and what to do about it (ex. they need to cry to let us know) | ||
| This is another one I would use the word "primal" to describe | ||
| Watching other's reaction to painful situations | ||
| Whether or not they could ignore it or recall feeling pain before and cope | ||
| getting hungry | All children are capable of being hungry | |
| Biological Need | ||
| child is controlling himself | ||
| everyone gets hungry no matter what the age is . it depends on how we show we are hungry | ||
| evolution, survival mechanisms | ||
| Getting hungry is one of the first things a newborn will experience | ||
| I think we're all just born with the feeling of hunger because he have to have nourishment | ||
| It's just survival at this point | ||
| It's simply built in | ||
| just something theyre born with | ||
| No, I think it's natural for all living beings to get hungry | ||
| None come to mind | ||
| the child being active and using enagy | ||
| The child has emotions that cause them to increase their desire to do this | ||
| the child just gets hungry- human nature | ||
| the more experience they have with food, the bigger the capacity for hunger | ||
| We are just preprogrammed to know when we are hungry. | ||
| Yeah, I think this one could be covered with just one answer and needs no further thought or explanation | ||
| You're born with this ability. Technically, you have this ability BEFORE you're born | ||
| Negative emotions | feeling distressed | again this is similar to the first one as they age they become more independent and spend less time distressed |
| All children can feel distressed at any given age | ||
| Babies show that they can be distressed from birth, or even during labor. | ||
| Distress is a biological survival mechanism | ||
| Distress is easily felt regardless of age. The factors change. A newborn gets distressed if they are hungry, or tired, or sick, etc. A 5 year old can get distressed if their favorite toy breaks | ||
| Feeling distressed is natural as the child has physical needs and cannot attend to them on their own | ||
| General living | ||
| I believe the feeling of distress is another innate ability that we are able to feel from the moment we're born. I think living, breathing creature that has feelings is the same way from birht | ||
| I do not imagine that any child of any age does not have this ability | ||
| I feel like feeling distressed can be attributed to human's evolution. Distress would be hard coded in the DNA (panic when something bad happens help survival) | ||
| I think babies crying shows distress | ||
| I think distress is an innate thing that's already there long before birth | ||
| I think it's similar to the pain feeling. They can feel distressed, but they may not understand the feeling | ||
| I think that feeling distressed is somewhat primal. Flight or flight, so to speak. I don't know if I would use the word "biological" to describe it, given my choice, but it's the closest out of the option provided. Then there are layers of socialization, etc. beyond that which amplify it | ||
| i think they always know distress in my opinion | ||
| if the child is left alone or isnt care for | ||
| It's an automatic response | ||
| It's the result of biological function, like pain response | ||
| Maybe hereditary? | ||
| Maybe the amount of nutrition a child gets throughout their life cycle | ||
| None come to mind | ||
| There could be physiological elements that prevent a child from being able to feel distress as acutely as other children | ||
| They are alive | ||
| They can feel stress from the people around them especially parents even inside the womb | ||
| You're born with all emotions, even if you don't understand them yet | ||
| feeling helpless | Again, you're born with emotions. You might not fully comprehend them, but you have the ability to feel them, nontheless | |
| as they age there is a tendency to become more independent and at that age less aware of danger | ||
| Feeling helpless is something that is learned through experience, not an innate ability | ||
| Feeling helplessness is not really an ability at all, it's an emotion. Most people are capable of feeling helpless | ||
| I don't think you can feel helpless until you become completely self aware, and I'm not sure when that happens | ||
| I feel like a kid would need to have some understanding of their emotions | ||
| I feel like this is something that a child knows | ||
| I feel that it changes. A younger child may not know or understand the feeling of helpless. An older child my feel they are invinsible or independent. | ||
| I really don't know what develops a child's sense to feel helpless | ||
| I think helplessness stems from a fundamental fear of losing something we've got or not getting something we want (regardless of how well we are able to articulate or understand what that is in a sophisticated way) | ||
| I think it is just mental development. Young children don't really feel self pity | ||
| I think that it all surrounds object perception/realization. Once the child realizes that "food" makes them not hungry or "moving" makes them not hurt | ||
| I think this is something that happens later in life | ||
| My zero is because at the beginning, there's no feeling of self, so you can't feel helpless if you don't feel the self (in my opinion) | ||
| None come to mind | ||
| Not being taken care of or given what it wants would eventually make them real familiar with feeling helpless | ||
| The child has negative experiences that increase this | ||
| The child notices when it is not cared for properly by feelings of hunger,pain etc | ||
| The child's need to be loved and cared for contributes to their feeling helpless, since they are generally moreso helpless in the earlier stages of life than later | ||
| The child's needs change | ||
| This is also hard to rate but I go with my same rationale for happiness. Infants ARE helpless and I think on some leel they know it | ||
| Very complicated question, as feeling helpless first requires a concept of helplessness (and in a way a name for that condition.) It might require the knowledge of help itself. So if a creature has never known help of any kind how can they know helplessness. Also there is learned helplessness which is different from experiencing the feeling of being helpless. I almost need more specific info to answer this one. But I just answered in terms of general capability | ||
| Social abilities & positive emotions | feeling happy | A child that feels love feels happiness |
| Again, I think happiness (albeit, not amusement) is something relatively innate related to security & safety. I feel like this should be an incredibly easy question to pinpoint to a specific cause, and yet somehow it is not at all | ||
| As the brain grows, we are capable of sensating new emotions. | ||
| Brain development in general plays a key role, as the child gets older, the brain is able to comprehend emotion | ||
| Children can feel joy at a certain age | ||
| I suppose in a way opportunities or reasons to be happy play a role in the development of capacity for happiness | ||
| I think it's innate in children..Sometimes they interpret it differently as they get older though | ||
| I think its just within them to feel happy | ||
| I think the ability to know they are happy also goes along with brain development | ||
| it is a normal emotion | ||
| None come to mind | ||
| Other external stimuli encourages this | ||
| social environment | ||
| The child becomes more able to notice their own feelings | ||
| the child senses the love shared around him\her | ||
| The child's interests/needs change | ||
| This is a harder one for me to rate. I would say some form of emotional range is present at birth. It would have to be. Beyond that I do not know | ||
| Well fed or not | ||
| yes it depends if the child is in a good home or a bad home | ||
| You're born with this ability | ||
| learning from other people | after the age of 2 or 3 the ability to learn langauges begins to diminish hence my ratings after that age | |
| Children have the capacity to learn from people at a certain age | ||
| I feel like this is just the nature of human beings to learn from the others around (we wouldn't be where we are without it) | ||
| I think there are a ton of factors that lend to learning. I think it was one of the major purposes of the brain and so many, many parts of the brain and body contribute to it | ||
| kids love to copy people they are around and that how they learn | ||
| Learning is an innate ability. Generally, everyone is born with the ability to learn | ||
| None come to mind | ||
| Relationship with parents | ||
| The amount of nutrition the child gets throughout their life cycle | ||
| The more they age, the more connections the brain makes. | ||
| Cognition & control | controlling their emotions | A 0 to 9 month old infant is never wrong. You don't tell a 0 to 9 month old infant "No" or "Not right now". After 9 months, then you can start having that conversation |
| Children can learn to control their emotions after a certain age | ||
| Children imitate other children | ||
| Environmental influences | ||
| having very good role models and a good home to learn these in | ||
| I don't really think there is anything else. It is taught and seen | ||
| I don't think anyone can control thier emotions (make yourself sad, NOW!), only the expression of them | ||
| I know a lot of adults that can't control their emotions. This is something you either have or don't have, | ||
| I would say the type of parenting being done by the parents | ||
| I'm afraid I don't have much to add here. Most of your questions covered it. I would say their environment, but that was, more or less, covered in the questions | ||
| Learning from adults reaction to behavior | ||
| Learning that outcomes can change depending on the the emotion they display | ||
| Maybe an innate connection with being, life, presence, consciousness (that thing which can't quite be named.) | ||
| None come to mind | ||
| overall social environment | ||
| The child has life experiences after birth that increase their ability to do this | ||
| They aren't necessarily learning to control the emotion so much as the expression | ||
| They brains haven't fully developed they cab;'t help themselves and are still exploring their emotions | ||
| They react to positive or negative responses from other humans, if they control, or don't control their emotions | ||
| watching tv or movies | ||
| reasoning about things | Again, up to about 9 months old the wiring is just not there. It beginbs to develop but... to be honest reasoning is something many adults fail at | |
| Genetic driven intelligence ability | ||
| injury/disability | ||
| It's mostly just natural mental development | ||
| Maybe the amount of food/nutrition they get throughout their life cycle | ||
| No not that I'm aware of | ||
| None come to mind | ||
| not sure if they can have those skills | ||
| Nuerological development can impact this..for example if they develop brain injuries...autism etc | ||
| Reasoning I feel is something they learn as the brain develops | ||
| The more they age, the more they comprehend and understand. | ||
| watching tv or videos |
Plots: Most important developmental factor, by domain and capacity
ggplot(d_dev_factor_most_important_choice %>%
mutate(domain = recode_factor(
domain,
"BOD" = "Bodily sensations",
"NEG" = "Negative emotions",
"POS" = "Social abilities & positive emotions",
"COG" = "Cognition & control",
.default = NA_character_)) %>%
mutate(dev_factor = case_when(
grepl("Other", response) ~ "other",
grepl("teach", response) ~ "people teach",
grepl("experiments", response) ~ "experiments",
grepl("womb", response) ~ "womb experiences",
grepl("interacts", response) ~ "interacts people",
grepl("preprogrammed", response) ~ "preprogrammed",
grepl("objects", response) ~ "observes objects",
grepl("observes the people", response) ~ "observes people",
grepl("body grows", response) ~ "body grows",
grepl("brain changes", response) ~ "brain changes",
grepl("senses improve", response) ~ "senses improve"),
dev_factor = factor(
dev_factor,
levels = c(gsub("_", " ",
levels(d_dev_factor_rating$dev_factor)),
"other"))) %>%
mutate_at(vars(capacity, response),
funs(gsub("_", " ", .))),
aes(x = dev_factor, fill = domain)) +
facet_wrap(~ domain ~ capacity, ncol = 4) +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
geom_bar() +
labs(title = "Choice of most important developmental factor (Version 1)",
x = "Developmental factor",
y = "Number of participants")
d_dev_factor_most_important_choice_other <- d_dev_factor_most_important_choice %>%
filter(grepl("Other", response)) %>%
count(domain, capacity)
ggplot(d_dev_factor_most_important_choice %>%
mutate(domain = recode_factor(
domain,
"BOD" = "Bodily sensations",
"NEG" = "Negative emotions",
"POS" = "Social abilities & positive emotions",
"COG" = "Cognition & control",
.default = NA_character_)) %>%
mutate(dev_factor = case_when(
grepl("Other", response) ~ "other",
grepl("teach", response) ~ "people teach",
grepl("experiments", response) ~ "experiments",
grepl("womb", response) ~ "womb experiences",
grepl("interacts", response) ~ "interacts people",
grepl("preprogrammed", response) ~ "preprogrammed",
grepl("objects", response) ~ "observes objects",
grepl("observes the people", response) ~ "observes people",
grepl("body grows", response) ~ "body grows",
grepl("brain changes", response) ~ "brain changes",
grepl("senses improve", response) ~ "senses improve"),
dev_factor = factor(
dev_factor,
levels = c(gsub("_", " ",
levels(d_dev_factor_rating$dev_factor)),
"other"))) %>%
mutate_at(vars(capacity, response),
funs(gsub("_", " ", .))) %>%
filter(dev_factor != "other"),
aes(x = reorder(capacity, as.numeric(domain)), fill = domain)) +
facet_wrap(~ dev_factor, ncol = 5) +
theme(legend.position = "top",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
geom_bar() +
labs(title = "Choice of most important developmental factor (Version 2)",
subtitle = paste0("Excluding answers of 'other' (n < ",
max(d_dev_factor_most_important_choice_other$n) + 1, ") for any capacity"),
color = "Domain",
x = "Capacity (by domain)",
y = "Number of participants")
d_dev_factor_most_important_free %>%
mutate(domain = recode_factor(
domain,
"BOD" = "Bodily sensations",
"NEG" = "Negative emotions",
"POS" = "Social abilities & positive emotions",
"COG" = "Cognition & control",
.default = NA_character_),
capacity = gsub("_", " ", capacity)) %>%
distinct(domain, capacity, response) %>%
arrange(domain, capacity, response) %>%
kable(caption = "Free responses when 'other' was selected as most important developmental factor") %>%
kable_styling() %>%
collapse_rows(1:2)| domain | capacity | response |
|---|---|---|
| Bodily sensations | feeling pain | everyone is capable of feeling pain. |
| getting hungry | All children are capable of being hungry. | |
| Negative emotions | feeling distressed | As soon as his consious develops they are considered alive and capable of feeling. |
| emotional development | ||
| most id not all children are capable of experinceing distress not matter what. | ||
| feeling helpless | Experiences in life | |
| I really don't know. | ||
| Not sure helplessness is felt so young | ||
| Social abilities & positive emotions | feeling happy | Every child is capable of feeling happiness no matter what. |
| it is just a normal response too something enjoyed | ||
| learning from other people | not sure | |
| Cognition & control | controlling their emotions | impossibe |
Relationships between developmental factor ratings and capacity ratings
temp <- d_cap_rating %>%
rename(response_cap_rating = response) %>%
full_join(d_dev_factor_rating_byclust %>%
rename(response_dev_factor = response))Joining, by = c("ResponseId", "domain", "capacity")ggplot(temp %>% filter(target == "target00mo"),
aes(x = response_dev_factor, y = response_cap_rating,
color = domain, group = capacity)) +
facet_grid(cols = vars(domain, capacity),
rows = vars(dev_factor_cluster, dev_factor)) +
geom_jitter(alpha = 0.2) +
geom_smooth(color = "black", method = "lm") +
theme(legend.position = "none") +
labs(title = "Relationship between endorsement of extrinsic vs. intrinsic developmental factors and attributions of capaciites at birth",
x = "Endorsement of developmental factor",
y = "Capacity attriubtion")
ggplot(temp,
aes(x = response_dev_factor, y = response_cap_rating)) +
facet_grid(dev_factor_cluster ~ target_ord) +
geom_jitter(alpha = 0.01, color = "gray") +
geom_smooth(aes(color = dev_factor, fill = dev_factor), method = "lm") +
scale_color_brewer(palette = "Paired") +
scale_fill_brewer(palette = "Paired") +
theme(legend.position = "bottom") +
labs(title = "Relationship between endorsement of extrinsic vs. intrinsic developmental factors and attributions of capacities at various target ages",
x = "Endorsement of developmental factor",
y = "Capacity attriubtion")
ggplot(temp %>% filter(domain == "BOD"),
aes(x = response_dev_factor, y = response_cap_rating)) +
facet_grid(rows = vars(capacity, dev_factor_cluster),
cols = vars(target_ord)) +
geom_jitter(alpha = 0.04, color = "gray") +
geom_smooth(aes(color = dev_factor, fill = dev_factor), method = "lm") +
scale_color_brewer(palette = "Paired") +
scale_fill_brewer(palette = "Paired") +
theme(legend.position = "bottom") +
labs(title = "BODILY SENSATIONS: Relationship between endorsement of extrinsic vs. intrinsic developmental factors and attributions of capacities at various target ages",
x = "Endorsement of developmental factor",
y = "Capacity attriubtion")
ggplot(temp %>% filter(domain == "NEG"),
aes(x = response_dev_factor, y = response_cap_rating)) +
facet_grid(rows = vars(capacity, dev_factor_cluster),
cols = vars(target_ord)) +
geom_jitter(alpha = 0.04, color = "gray") +
geom_smooth(aes(color = dev_factor, fill = dev_factor), method = "lm") +
scale_color_brewer(palette = "Paired") +
scale_fill_brewer(palette = "Paired") +
theme(legend.position = "bottom") +
labs(title = "NEGATIVE EMOTIONS: Relationship between endorsement of extrinsic vs. intrinsic developmental factors and attributions of capacities at various target ages",
x = "Endorsement of developmental factor",
y = "Capacity attriubtion")
ggplot(temp %>% filter(domain == "POS"),
aes(x = response_dev_factor, y = response_cap_rating)) +
facet_grid(rows = vars(capacity, dev_factor_cluster),
cols = vars(target_ord)) +
geom_jitter(alpha = 0.04, color = "gray") +
geom_smooth(aes(color = dev_factor, fill = dev_factor), method = "lm") +
scale_color_brewer(palette = "Paired") +
scale_fill_brewer(palette = "Paired") +
theme(legend.position = "bottom") +
labs(title = "POSITIVE EMOTIONS & SOCIAL ABILITIES: Relationship between endorsement of extrinsic vs. intrinsic developmental factors and attributions of capacities at various target ages",
x = "Endorsement of developmental factor",
y = "Capacity attriubtion")
Social abilities & positive emotions
Linear effects only
Non-linear effects